#ifndef __VECTOR3D_H__
#define __VECTOR3D_H__
#include <assert.h>
#include <math.h>
#include <float.h>
using namespace std;
class Vector3D{
        private:
        double v[3];
        public:
                Vector3D();
                Vector3D(double x, double y, double z);
                Vector3D(const Vector3D& o);
                ~Vector3D();
                Vector3D& operator=(const Vector3D& o);
                Vector3D& operator+=(const Vector3D& o);
                Vector3D& operator-=(const Vector3D& o);
                const Vector3D& operator-() const;
                const Vector3D& operator+(const Vector3D& o) const;
                const Vector3D& operator-(const Vector3D& o) const;
                const Vector3D operator*(const double f) const;
                const Vector3D operator/(const double f) const;
                double& operator[](int i);
                const double operator*(const Vector3D& o) const;
                const Vector3D crossProduct(const Vector3D& o) const;
                bool operator==(const Vector3D& o) const;
                bool operator!=(const Vector3D& o) const;
                const double norm() const;
                const double squareNorm() const;
                void normalize();
                friend ostream& operator<<(ostream& s, Vector3D& p);
                friend istream& operator>>(istream& s, Vector3D& p);
                friend const Vector3D operator*(double f, Vector3D p);
                friend const Vector3D operator/(double f, Vector3D p);
};

Vector3D::Vector3D() {
        v[0]=0;
        v[1]=0;
        v[2]=0;
}
Vector3D::Vector3D(double x, double y, double z){
        v[0] = x;
        v[1] = y;
        v[2] = z;
}
Vector3D::Vector3D(const Vector3D& o){
        for(int i=0;i<3;i++) v[i] = o.v[i];
}
Vector3D::~Vector3D(){ }
Vector3D& Vector3D::operator=(const Vector3D& o){
        if(&o != this)
                for(int i=0;i<3;i++) v[i] = o.v[i];
        return *this;
}
Vector3D& Vector3D::operator+=(const Vector3D& o){
        v[0]+=o.v[0];
        v[1]+=o.v[1];
        v[2]+=o.v[2];
        return *this;
}
Vector3D& Vector3D::operator-=(const Vector3D& o){
        v[0]-=o.v[0];
        v[1]-=o.v[1];
        v[2]-=o.v[2];
        return *this;
}
const Vector3D& Vector3D::operator-() const{
        return (Vector3D(0,0,0) - *this);
}
const Vector3D& Vector3D::operator+(const Vector3D& o) const{
        return Vector3D(*this)+=o;
}
const Vector3D& Vector3D::operator-(const Vector3D& o) const{
        return Vector3D(*this)-=o;
}
const Vector3D Vector3D::operator*(const double f) const {
        Vector3D temp = *this;
        temp.v[0] *=f;
        temp.v[1] *=f;
        temp.v[2] *=f;
        return temp;
}
const Vector3D Vector3D::operator/(const double f) const {
        Vector3D temp = *this;
        if(f==0) return Vector3D(0,0,0);
        double b = 1.0f / f;
        temp.v[0] *=b;
        temp.v[1] *=b;
        temp.v[2] *=b;
        return temp;
}
const double Vector3D::operator*(const Vector3D& o) const{
        return v[0]*o.v[0]+v[1]*o.v[1]+v[2]*o.v[2];
}
const inline Vector3D Vector3D::crossProduct(const Vector3D& o) const{
        return Vector3D(v[1]*o.v[2]-v[2]*o.v[1],v[0]*o.v[2]-v[2]*o.v[0],v[0]*o.v[1]-v[1]*o.v[0]);
}
bool Vector3D::operator==(const Vector3D& o) const{
        for(int i=0;i<3;i++) 
                if(fabs(v[i] - o.v[i]) > DBL_EPSILON) return false;     //Takes care of doubleing precision.
        return true;
}
bool Vector3D::operator!=(const Vector3D& o) const{
        return !((*this) == o);
}
const double Vector3D::norm() const{
        return sqrt(squareNorm());
}
const double Vector3D::squareNorm() const{
        return (*this)*(*this);
}
void Vector3D::normalize(){
        double squNorm = squareNorm();
        if(squNorm!=0){
                double _norm = sqrt(squNorm);
                v[0] /= _norm; v[1] /= _norm; v[2] /= _norm;
        }
}
ostream& operator<<(ostream& s, Vector3D& p){
        s<<p.v[0]<<" "<<p.v[1]<<" "<<p.v[2];
        return s;
}
istream& operator>>(istream& s, Vector3D& p){
        s>>p.v[0]>>p.v[1]>>p.v[2];
        return s;
}
const Vector3D operator*(const double f, const Vector3D p){
        return p*f;
}
const Vector3D operator/(const double f, const Vector3D p){
        return p/f;
}
double& Vector3D::operator[](int i){
        assert(i>-1 && i<3);
        return v[i];
}
#endif //__VECTOR3D_H__
